﻿#' @TODO 生成单因素筛选预后基因 森林图展示
#' @param unicox_res 单因素cox结果，需要包含gene HR p.value Hazard_Ratio lower_.95 upper_.95等表头
#' 
#' > head(univcox_res,5)
#'          gene       coef    p.value Hazard_Ratio lower_.95 upper_.95
#' SPRY2   SPRY2 -0.2238623 0.03005903    0.7994252 0.6530374 0.9786280
#' CCDC69 CCDC69 -0.2113007 0.02519092    0.8095306 0.6727964 0.9740536
#' PLK1     PLK1  0.1558183 0.04529280    1.1686139 1.0032713 1.3612054
#' TUBA1C TUBA1C  0.3104911 0.01054146    1.3640949 1.0752437 1.7305425
#' EZR       EZR  0.2933651 0.02835023    1.3409323 1.0315917 1.7430340
#'        logrank_pvalue wald_pvalue Likelihood_pvalue                  HR
#' SPRY2      0.03007443  0.03005903        0.02668130 0.799 (0.653-0.979)
#' CCDC69     0.02501043  0.02519092        0.02353125  0.81 (0.673-0.974)
#' PLK1       0.04497782  0.04529280        0.04598493 1.169 (1.003-1.361)
#' TUBA1C     0.01052236  0.01054146        0.01138811 1.364 (1.075-1.731)
#' EZR        0.02847897  0.02835023        0.02885348 1.341 (1.032-1.743)
#' @return 函数 forestplot 绘制图像结果
#' @param 
#' 
#' @Author WYK
#' 
univcox_geneforest <- function(univcox_res = NULL, output_dir = "./") {
  if (!"HR" %in% colnames(univcox_res)) {
    univcox_res$HR <- paste0(
      round(univcox_res$Hazard_Ratio, 3),
      " (",
      sprintf("%.3f", univcox_res$lower_.95),
      "-",
      sprintf("%.3f", univcox_res$upper_.95),
      ")"
    )
  }

  if (!dir.exists(sprintf("%soutput/univcox_gene_forest", output_dir))) {
    dir.create(sprintf("%soutput/univcox_gene_forest", output_dir), recursive = T)
  }

  tabletext <- cbind(
    c("Prognostic gene", univcox_res$gene),
    c("Hazard ratio (95% CI)", univcox_res$HR),
    c("P value", format.pval(univcox_res$p.value,digits = 1,eps = 0.001))
  )

  library(forestplot)

  min_num <- round(min(univcox_res$Hazard_Ratio) - .3, 1)
  max_num <- round(max(univcox_res$Hazard_Ratio) + 0.2, 1)

  tiff(
    filename = sprintf("%soutput/univcox_gene_forest/univcox_geneforest.tiff", output_dir),
    units = "px", width = 7.4 * 300, height = ncol(univcox_res) * .6 * 300, res = 300
  )
  forestplot(
    labeltext = tabletext, # 主要是以矩阵或者list形式将数据导入函数，最好以矩阵，因为数据一般都是矩阵的。
    mean = c(NA, univcox_res$Hazard_Ratio), # 误差条的均值
    lower = c(NA, univcox_res$lower_.95), # 误差条 95%置信区间下限
    upper = c(NA, univcox_res$upper_.95), # 误差条 95%置信区间上限
    clip = c(min_num, max_num),
    xlog = TRUE,
    zero = 1,
    boxsize = 0.35,
    hrzl_lines = list("2" = gpar(lwd = 2, col = "#99999922")),
    graphwidth = unit(70, "mm"),
    lineheight = "auto", # 行的高度，可以是数字，也可以是 unit 的形式
    # line.margin = .1, # 行与行之间的间隙的宽度
    colgap = unit(4, "mm"),
    lwd.zero = 2,
    lwd.ci = 2,
    cex = 0.9,
    ci.vertices = TRUE,
    ci.vertices.height = 0.2,
    col = fpColors(
      box = "#1c61b6",
      line = "#1c61b6",
      zero = "gray50"
    ), # 森林图横线以及点的颜色。
    # box：box（点估计值）的颜色
    # line：穿过方块的横线的颜色
    # zero：中间那条基准线的颜色
    # summary：summary中菱形的颜色
    # hrz_lines：表中第一条横线的颜色
    # eg：col=fpcolors(box=’royblue’,line=’darkblue’, summary=’royblue’, hrz_lines=’red’)
    xlab = "Hazard Ratio",
    graph.pos = 2, # 定位森林图所在的位置。通过数字来确定为第几列
    txt_gp = fpTxtGp(
      label = gpar(cex = 1.05),
      ticks = gpar(cex = 1),
      xlab = gpar(cex = 1.1),
    ), # 设置表格中文本的格式：用gpar进行赋值，其中cex为文本字体大小，ticks为坐标轴大小，xlab为坐标轴文字字体大小。
    # label：表格主体文字的格式
    # ticks：森林图下方的坐标轴的刻度文字格式
    # xlab：定义的x轴标题格式
    # title：标题文字的格式
    # eg：txt_gp=fpTxtGp(label=gpar(cex=1.25), ticks=gpar(cex=1.1), xlab=gpar(cex = 1.2), title=gpar(cex = 1.2))
    #
    # align = 'c', # 每列文字的对齐方式，偶尔会用到。如：align=c("l","c","c")l：左对齐r：右对齐c：居中对齐
    new_page = TRUE,
    xticks = c(min_num, 1, max_num)
  )
  dev.off()

  tiff(
    filename = sprintf("%soutput/univcox_gene_forest/univcox_geneforest_dpi72.tiff", output_dir),
    units = "px", width = 7.4 * 90, height = ncol(univcox_res) * .6 * 90, res = 90
  )
  forestplot(
    labeltext = tabletext,
    mean = c(NA, univcox_res$Hazard_Ratio),
    lower = c(NA, univcox_res$lower_.95),
    upper = c(NA, univcox_res$upper_.95),
    clip = c(round(min(univcox_res$Hazard_Ratio) + 0.2, 1), round(max(univcox_res$Hazard_Ratio) + 0.2, 1)),
    xlog = TRUE,
    zero = 1,
    boxsize = 0.35,
    hrzl_lines = list("2" = gpar(lwd = 2, col = "#99999922")),
    graphwidth = unit(70, "mm"),
    lineheight = "auto",
    colgap = unit(4, "mm"),
    lwd.zero = 2,
    lwd.ci = 2,
    cex = 0.9,
    ci.vertices = TRUE, ci.vertices.height = 0.2,
    col = fpColors(
      box = "#1c61b6",
      line = "#1c61b6",
      zero = "gray50"
    ),
    xlab = "Hazard Ratio",
    graph.pos = 2,
    txt_gp = fpTxtGp(
      label = gpar(cex = 1.05),
      ticks = gpar(cex = 1),
      xlab = gpar(cex = 1.1),
    ),
    new_page = TRUE,
    xticks = c(min_num, 1, max_num)
  )
  dev.off()

  pdf(file = sprintf("%soutput/univcox_gene_forest/univcox_geneforest.pdf", output_dir), width = 7.4, height = ncol(univcox_res) * .6)
  forestplot(
    labeltext = tabletext,
    mean = c(NA, univcox_res$Hazard_Ratio),
    lower = c(NA, univcox_res$lower_.95),
    upper = c(NA, univcox_res$upper_.95),
    clip = c(round(min(univcox_res$Hazard_Ratio) + 0.2, 1), round(max(univcox_res$Hazard_Ratio) + 0.2, 1)),
    xlog = TRUE,
    zero = 1,
    boxsize = 0.25,
    hrzl_lines = list("2" = gpar(lwd = 2, col = "#99999922")),
    graphwidth = unit(70, "mm"),
    lineheight = "auto",
    colgap = unit(4, "mm"),
    lwd.zero = 2,
    lwd.ci = 2,
    cex = 0.9,
    ci.vertices = TRUE,
    ci.vertices.height = 0.2,
    col = fpColors(
      box = "#1c61b6",
      line = "#1c61b6",
      zero = "gray50"
    ),
    xlab = "Hazard Ratio",
    graph.pos = 2,
    txt_gp = fpTxtGp(
      label = gpar(cex = 1.05),
      ticks = gpar(cex = 1),
      xlab = gpar(cex = 1.1),
    ),
    new_page = F,
    xticks = c(min_num, 1, max_num)
  )
  dev.off()
}

